Marginalizing shocks
Cover that shocks which I ought not to see.
1 Model reparameterization
We marginalize the shocks in the mixture model! Going from:
# Non-shocky
target += normal_lpdf(y[n] | mu[n] + delta_sck[n], sigma)
target += normal_lpdf(delta_sck[n] | 0, inner_tau_sck)
# Shocky
target += normal_lpdf(y[n] | mu[n] + delta_sck[n], sigma)
target += normal_lpdf(delta_sck[n] | 0, outer_tau_sck)to a log-likelihood where shocks are marginalized out:
# Non-shocky
target += normal_lpdf(y[n] | mu[n], eta)
# Shocky
target += normal_lpdf( y[n] | mu[n], eta * sqrt{1 + nu / eta^2})where: \eta = \sqrt{ \sigma^2 + {\tau_\text{inner}^\text{shock}}^2 } \nu = {\tau_\text{outer}^\text{shock}}^2-{\tau_\text{inner}^\text{shock}}^2 In Stan, we now have in the model block:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks.stan'))[230:315])model {
alpha ~ normal(0, 0.69); // 0.2 mm <~ exp(alpha) * 1 mm <~ 5 mm
beta_gdd ~ normal(0, log(1.8) / 2.57); // -log(1.8) <~ beta_gsl <~ log(1.8)
beta_sm ~ normal(0, log(1.8) / 2.57); // -log(1.8) <~ beta_gsl <~ log(1.8)
beta_vpd ~ normal(0, log(1.8) / 2.57); // -log(1.8) <~ beta_gsl <~ log(1.8)
rho_sp ~ lognormal(2.65, 0.135); // 10 <~ rho <~ 20
gamma_sp ~ normal(0, log(10) / 2.57); // 0 <~ gamma <~ log(10)
for (s in 1:N_stands)
f_tilde_sh[s] ~ normal(0, 1);
rho_sh ~ lognormal(1.7, 0.26); // 3 <~ rho_sh <~ 10
gamma_sh ~ normal(0, log(3) / 2.57); // 0 <~ gamma_sh <~ log(3)
eta ~ gamma(6, 43.5);
nu ~ gamma(0.5, 0.033);
phi_sck ~ beta(2, 20);
omega_conc_sck ~ beta(230, 14); // 0.9 <~ omega_conc_sck <~ 0.97 (most trees, but not ALL trees)
omega_nonconc_sck ~ beta(1, 20); // 0 <~ omega_conc_sck <~ 0.15 (should be rare, but... who knows?)
vector[N] mu;
for (t in 1:N_trees) {
array[N_years[t]] int tree_idxs
= linspaced_int_array(N_years[t],
tree_start_idxs[t],
tree_end_idxs[t]);
array[N_years[t]] int all_years_idxs_tree = all_years_idxs[tree_idxs];
array[N_years[t]] real years_tree = years[tree_idxs];
int stand_idx = stand_idxs[t];
int species_idx = species_idxs[t];
vector[N_years[t]] gdd_obs_tree = gdd_obs[tree_idxs];
vector[N_years[t]] sm_obs_tree = sm_obs[tree_idxs];
vector[N_years[t]] vpd_obs_tree = vpd_obs[tree_idxs];
mu[tree_idxs] = alpha
+ beta_gdd[species_idx] * (gdd_obs_tree - gdd0)
+ beta_sm[species_idx] * (sm_obs_tree - sm0)
+ beta_vpd[species_idx] * (vpd_obs_tree - vpd0)
+ f_sh[stand_idx, all_years_idxs_tree];
f_tilde[tree_idxs] ~ normal(0,1);
}
// Observational model with marginalized shocks
for (s in 1:N_stands) {
array[N_stand_trees[s]] int stand_trees_idxs = linspaced_int_array(N_stand_trees[s],
stand_trees_start_idxs[s], stand_trees_end_idxs[s]);
vector[N_stand_years[s]] log_p0 = rep_vector(0, N_stand_years[s]);
vector[N_stand_years[s]] log_p1 = rep_vector(0, N_stand_years[s]);
for(ts in 1:N_stand_trees[s]){
int t = stand_trees_idxs[ts];
array[N_years[t]] int tree_idxs = linspaced_int_array(N_years[t], tree_start_idxs[t], tree_end_idxs[t]);
array[N_years[t]] int all_years_idxs_tree = all_years_idxs[tree_idxs];
for(y in 1:N_years[t]) {
int ys = all_years_idxs_tree[y]-stand_start_years_idxs[s]+1;
log_p0[ys] += log_mix(omega_nonconc_sck,
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]]
+ f[tree_idxs[y]], eta * sqrt(1 + nu / eta^2)),
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]]
+ f[tree_idxs[y]], eta));
log_p1[ys] += log_mix(omega_conc_sck,
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]]
+ f[tree_idxs[y]], eta * sqrt(1 + nu / eta^2)),
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]]
+ f[tree_idxs[y]], eta));
}
}
for(y in 1:N_stand_years[s]) {
target += log_mix(phi_sck[s], log_p1[y], log_p0[y]);
}
}
}2 But… my shocks…
Our large shocks are likely tied to interesting ecological stuff… But by marginalizing them, we don’t have inferences on them anymore.
It could be doable to “to marginalize out \delta_\text{shock} in the “small shock” branch of the mixture model but leave \delta_\text{shock} explicit in the “large shock” branch of the mixture model” (Mike).
Since I’m not totally sure how to do that, I decided to sample the \delta_\text{shock} afterwards, to recover their posterior distributions in the generated quantities block. This is an ongoing project, so for now I recover only the shock state (z_t).
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks.stan'))[351:405])generated quantities {
array[N] int sck_state; // latent state, zt = 0 or zt = 1
array[N] real log_rw_pred;
for (t in 1:N_trees) {
array[N_years[t]] int tree_idxs
= linspaced_int_array(N_years[t],
tree_start_idxs[t],
tree_end_idxs[t]);
array[N_years[t]] int all_years_idxs_tree = all_years_idxs[tree_idxs];
array[N_years[t]] real years_tree = years[tree_idxs];
int stand_idx = stand_idxs[t];
int species_idx = species_idxs[t];
vector[N_years[t]] gdd_obs_tree = gdd_obs[tree_idxs];
vector[N_years[t]] sm_obs_tree = sm_obs[tree_idxs];
vector[N_years[t]] vpd_obs_tree = vpd_obs[tree_idxs];
vector[N_years[t]] mu = alpha
+ beta_gdd[species_idx] * (gdd_obs_tree - gdd0)
+ beta_sm[species_idx] * (sm_obs_tree - sm0)
+ beta_vpd[species_idx] * (vpd_obs_tree - vpd0)
+ f_sh[stand_idx, all_years_idxs_tree];
// mixture weight for shock
real mw_shock = phi_sck[stand_idx]*omega_conc_sck + (1-phi_sck[stand_idx])*omega_nonconc_sck;
for(y in 1:N_years[t]){
// mixture components
real log_pshock = log(mw_shock) + normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[y] + f[tree_idxs[y]], eta * sqrt(1 + nu / eta^2));
real log_pshock_plus_pnonshock = log_mix(mw_shock,
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[y] + f[tree_idxs[y]], eta * sqrt(1 + nu / eta^2)),
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[y] + f[tree_idxs[y]], eta));
// probability of shock state
real lambda_shock = exp(log_pshock - log_pshock_plus_pnonshock);
sck_state[tree_idxs[y]] = bernoulli_rng(lambda_shock); // or something like categorical_rng(lambda_shock);?
if(sck_state[tree_idxs[y]] == 1) {
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]], eta * sqrt(1 + nu / eta^2));
}else{
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]], eta);
}
}
}
}3 Test with a small dataset
Let’s run this new model on a very small dataset, 21 trees and 1 site.
The good thing here is that I don’t need to tune the (non-)centered parametrization anymore.
data <- readRDS(file = file.path(wd, 'output', 'model', 'data_15nov2025_az592.rds'))
data$N_stand_trees <- array(data$N_stand_trees, dim = 1)
data$N_stand_years <- array(data$N_stand_years, dim = 1)
data$stand_start_years_idxs <- array(data$stand_start_years_idxs, dim = 1)
data$stand_trees_end_idxs <- array(data$stand_trees_end_idxs, dim = 1)
data$stand_trees_start_idxs <- array(data$stand_trees_start_idxs, dim = 1)
if(FALSE){
fit <- stan(file=file.path(wd, 'model', 'stan', 'model9_marginalshocks.stan'),
data=data, seed=5838293, chains = 4,
warmup=1000, iter=2024, refresh=10)
saveRDS(fit, file.path(wd, 'model/shocks/output/marginal', 'fit.rds'))
}else{
fit <- readRDS(file.path(wd, 'model/shocks/output/marginal', 'fit.rds'))
}Let’s look at the diagnostics:
diagnostics <- util$extract_hmc_diagnostics(fit)
util$check_all_hmc_diagnostics(diagnostics) All Hamiltonian Monte Carlo diagnostics are consistent with reliable
Markov chain Monte Carlo.samples <- util$extract_expectand_vals(fit)
base_samples <- util$filter_expectands(samples,
c('alpha',
'beta_gdd', 'beta_sm', 'beta_vpd',
'rho_sh', 'gamma_sh',
'rho_sp', 'gamma_sp',
'eta', 'nu',
'phi_sck', 'omega_conc_sck', 'omega_nonconc_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)All expectands checked appear to be behaving well enough for reliable
Markov chain Monte Carlo estimation.Clean! Let’s continue.
These are the latent shock states:
set.seed(123456)
random_trees <- sample(1:data$N_trees, 4)
par(mfrow = c(2,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('sck_state[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=paste0("Shock state - tree ",t),
display_ylim=c(-0.1, 1.1), display_xlim = c(1980, 2010))
}
4 Next step
Is it possible to recover the posterior of \delta_\text{shock}(k) ~|~ y(k), without having {\tau_\text{inner}^\text{shock}} or {\tau_\text{outer}^\text{shock}}?
Maybe a first step could be to fix: {\tau_\text{inner}^\text{shock}} = 0?
Because for now my predictions look like:
par(mfrow = c(2,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=paste0("Shock state - tree ",t),
display_ylim=c(-10, 5), display_xlim = c(1980, 2010))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.7, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.4, col="black")
}
5 Assuming that {\tau_\text{inner}^\text{shock}} = 0
Note
Update Nov. 26 10pm: I tried another method in the section 8 to marginalize only the non-shocks, and reconstruct a posteriori \delta_\text{shock} that are \neq 0 only for large shocks. This does not use the normal-normal conjugancy anymore.
If I assume that {\tau_\text{inner}^\text{shock}} = 0, I can reconstruct the large shocks using the normal-normal conjugancy (merci Google): \delta_\text{shock}(k) ~|~ y_k \sim \text{normal} \left( ({\tau_\text{outer}^\text{shock}}^2 / ({\tau_\text{outer}^\text{shock}}^2 + \sigma^2)) \times y_k ~,~ \sqrt{({\tau_\text{outer}^\text{shock}}^2 \times \sigma^2) / ({\tau_\text{outer}^\text{shock}}^2 + \sigma^2)}\right)
In the generated quantities block, we can then have:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks_zerotauinner.stan'))[398:405]) if(sck_state[tree_idxs[y]] == 1) {
// we can reconstruct shock posterior using the normal-normal conjugancy
real conjugate_mean = (outer_tau_sck^2 / (outer_tau_sck^2 + sigma^2)) * log_rw_obs[tree_idxs[y]];
real conjugate_sd = sqrt((outer_tau_sck^2 * sigma^2) / (outer_tau_sck^2 + sigma^2));
delta_sck[tree_idxs[y]] = normal_rng(conjugate_mean, conjugate_sd);
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]] + delta_sck[tree_idxs[y]], sigma);
}else{
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]], sigma);if(FALSE){
fit2 <- stan(file=file.path(wd, 'model', 'stan', 'model9_marginalshocks_zerotauinner.stan'),
data=data, seed=5838293, chains = 4,
warmup=1000, iter=2024, refresh=10)
saveRDS(fit2, file.path(wd, 'model/shocks/output/marginal', 'fit2.rds'))
}else{
fit2 <- readRDS(file.path(wd, 'model/shocks/output/marginal', 'fit2.rds'))
}Diagnostics are cleaned:
diagnostics <- util$extract_hmc_diagnostics(fit2)
util$check_all_hmc_diagnostics(diagnostics) All Hamiltonian Monte Carlo diagnostics are consistent with reliable
Markov chain Monte Carlo.samples <- util$extract_expectand_vals(fit2)
base_samples <- util$filter_expectands(samples,
c('alpha',
'beta_gdd', 'beta_sm', 'beta_vpd',
'rho_sh', 'gamma_sh',
'rho_sp', 'gamma_sp',
'sigma', 'outer_tau_sck',
'phi_sck', 'omega_conc_sck', 'omega_nonconc_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)All expectands checked appear to be behaving well enough for reliable
Markov chain Monte Carlo estimation.We can look at the latent shocks:
par(mfrow = c(2,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('delta_sck[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=paste0("Delta shock - tree ",t),
display_ylim=c(-10, 0.1), display_xlim = c(1980, 2010))
}
And at the predictions:
par(mfrow = c(2,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=paste0("log(rw) - tree ",t),
display_ylim=c(-12, 1), display_xlim = c(1980, 2010))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.7, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.4, col="black")
}
However, it seems that the shocks are bit too large…
par(mfrow = c(2,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
baseline_values = data$log_rw_obs[idxs_t],
xlab="", ylab=paste0("Residuals - tree ",t),
display_ylim=c(-3, 1), display_xlim = c(1980, 2010),
residual = TRUE)
}
And, let’s not forget the stand-level GP, which is still wiggly:
par(mfrow = c(1,1), mar = c(2, 5, 3, 1))
s <- 1
names <- sapply(1:data$N_all_years,
function(y) paste0('f_sh[', s, ',', y, ']'))
util$plot_realizations(samples, names, plot_xs = data$all_years[1:data$N_all_years], N_plots = 50,
main = paste0('Stand ', data$uniq_stand_ids[s]), ylab = 'f_sh',
display_xlim = data$all_years[c(1,data$N_all_years)], display_ylim = c(-2,2))
6 More data!
Let’s expand this stuff to moooore data!
914 trees, this take more than 7 hours to run.
data <- readRDS(file = file.path(wd, 'output/model', 'data_15nov2025_azPIPO.rds'))
if(FALSE){
fit3 <- stan(file=file.path(wd, 'model', 'stan', 'model9_marginalshocks_zerotauinner.stan'),
data=data, seed=5838293, chains = 4,
warmup=1000, iter=2024, refresh=10)
saveRDS(fit3, file.path(wd, 'model/shocks/output/marginal', 'fit3.rds'))
}else{
fit3 <- readRDS(file.path(wd, 'model/shocks/output/marginal', 'fit3.rds'))
}diagnostics <- util$extract_hmc_diagnostics(fit3)
util$check_all_hmc_diagnostics(diagnostics) All Hamiltonian Monte Carlo diagnostics are consistent with reliable
Markov chain Monte Carlo.samples <- util$extract_expectand_vals(fit3)
base_samples <- util$filter_expectands(samples,
c('alpha',
'beta_gdd', 'beta_sm', 'beta_vpd',
'rho_sh', 'gamma_sh',
'rho_sp', 'gamma_sp',
'sigma', 'outer_tau_sck',
'phi_sck', 'omega_conc_sck', 'omega_nonconc_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)alpha:
Chain 2: hat{ESS} (94.815) is smaller than desired (100).
Chain 3: hat{ESS} (96.226) is smaller than desired (100).
beta_gdd[1]:
Chain 2: hat{ESS} (83.445) is smaller than desired (100).
beta_sm[1]:
Chain 2: hat{ESS} (70.945) is smaller than desired (100).
rho_sp[1]:
Chain 1: hat{ESS} (76.078) is smaller than desired (100).
Chain 2: hat{ESS} (37.400) is smaller than desired (100).
Chain 3: hat{ESS} (56.188) is smaller than desired (100).
Chain 4: hat{ESS} (65.237) is smaller than desired (100).
sigma:
Chain 2: hat{ESS} (38.755) is smaller than desired (100).
Chain 3: hat{ESS} (36.254) is smaller than desired (100).
Chain 4: hat{ESS} (89.635) is smaller than desired (100).
outer_tau_sck:
Chain 3: hat{ESS} (46.249) is smaller than desired (100).
phi_sck[12]:
Chain 3: hat{ESS} (68.030) is smaller than desired (100).
omega_conc_sck:
Chain 2: Right tail hat{xi} (0.871) exceeds 0.25.
Chain 4: Right tail hat{xi} (0.657) exceeds 0.25.
Split hat{R} (1.685) exceeds 1.1.
Chain 2: hat{ESS} (7.744) is smaller than desired (100).
Chain 3: hat{ESS} (4.191) is smaller than desired (100).
Chain 4: hat{ESS} (10.152) is smaller than desired (100).
omega_nonconc_sck:
Split hat{R} (1.353) exceeds 1.1.
Chain 2: hat{ESS} (10.463) is smaller than desired (100).
Chain 3: hat{ESS} (5.389) is smaller than desired (100).
Chain 4: hat{ESS} (15.782) is smaller than desired (100).
Large tail hat{xi}s suggest that the expectand might not be
sufficiently integrable.
Split Rhat larger than 1.1 suggests that at least one of the Markov
chains has not reached an equilibrium.
Small empirical effective sample sizes result in imprecise Markov chain
Monte Carlo estimators.The model clearly wants to concentrate on lower values of \omega_\text{shock}^\text{concordant}.
util$plot_pairs_by_chain(samples[['omega_conc_sck']], 'omega_conc_sck',
samples[['omega_nonconc_sck']], 'omega_nonconc_sck')
If we look at all shock-related parameters, we can see that both modes for \omega_\text{shock}^\text{concordant} are very far from our prior domain expertise.
# Quickly look at posteriors
nf <- layout(
matrix(c(1,1,1,2,2,2,
3,3,4,4,5,5),
ncol=6, byrow=TRUE)
)
util$plot_expectand_pushforward(samples[['sigma']], 50,
flim = c(0,0.5),
display_name=expression(tau[shock]^inner))
xs <- seq(0, 2, 0.01)
ys <- dnorm(xs, 0, log(1.1) / 2.57)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples[['outer_tau_sck']], 50,
flim = c(0,10),
display_name=expression(tau[shock]^outer))
xs <- seq(0, 10, 0.01)
ys <- dnorm(xs, 0, 10 / 2.57)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples[['phi_sck[1]']], 50,
flim = c(0,0.5),
display_name=expression(phi[schock]))
for(s in 2:data$N_stands){
util$plot_expectand_pushforward(samples[[paste0('phi_sck[',s,']')]], 50,
flim = c(0,0.5),
display_name=expression(phi[schock]), add = T)
}
xs <- seq(0, 0.5, 0.01)
ys <- dbeta(xs, 2, 20)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples[['omega_conc_sck']], 50,
flim = c(0,1),
display_name=expression(omega[shock]^concordant))
xs <- seq(0, 1, 0.01)
ys <- dbeta(xs, 230, 14)
lines(xs, ys, lwd=2.5, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples[['omega_nonconc_sck']], 50,
flim = c(0,0.5),
display_name=expression(omega[shock]^nonconcordant))
xs <- seq(0, 1, 0.01)
ys <- dbeta(xs, 1, 20)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
7 What are the concordant shocks captured by the model?
Let’s dig into these tree-level shocks (concordant, i.e. when the stand is shocking) that are less frequent that we initially thought.
We can first look at the stand we initially fitted alone (S14). Let’s look at 10 random trees.
trees <- which(grepl('az592' , data$uniq_tree_ids))
s1 <- unique(data$stand_idxs[trees])
random_trees <- sample(trees, 10)
par(mfrow = c(5,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=paste0("log(rw) - tree ",t),
display_ylim=c(-12, 1), display_xlim = c(1980, 2010))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.7, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.4, col="black")
}
mtext(paste0("Stand ", data$uniq_stand_ids[s1]), side = 3, line = -2, outer = TRUE, adj = 0)
On this stand, it does seem that most trees are shocking.
But on some other stand, the shock seem a bit less concordant across trees:
s2 <- 6
trees <- which(data$stand_idxs == s2)
random_trees <- sample(trees, 10)
par(mfrow = c(5,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=paste0("log(rw) - tree ",t),
display_ylim=c(-12, 1), display_xlim = c(1980, 2010))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.7, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.4, col="black")
}
mtext(paste0("Stand ", data$uniq_stand_ids[s2]), side = 3, line = -2, outer = TRUE, adj = 0)
The stand-level probability of shocks is not more important on the stand where the tree shocks are more concordant…
par(mfrow = c(1,1))
util$plot_expectand_pushforward(samples[[paste0('phi_sck[',s1,']')]],
50, flim = c(0,1), ylim = c(0,20),
display_name=expression(phi[shock]))
text(data$uniq_stand_ids[s1], x = 0.05, y = 16, col = util$c_dark)
util$plot_expectand_pushforward(samples[[paste0('phi_sck[',s2,']')]],
50, flim = c(0,1), ylim = c(0,20),
display_name=expression(phi[shock]), add = TRUE,
col = util$c_light_highlight)
text(data$uniq_stand_ids[s2], x = 0.17, y = 8, col = util$c_light_highlight)
Let’s look at another stand!
s3 <- 20
trees <- which(data$stand_idxs == s3)
random_trees <- sample(trees, 10)
par(mfrow = c(5,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=paste0("log(rw) - tree ",t),
display_ylim=c(-12, 1), display_xlim = c(1980, 2010))
abline(v = 2002, lty = 2, col = 'grey40')
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.7, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.4, col="black")
}
mtext(paste0("Stand ", data$uniq_stand_ids[s3]), side = 3, line = -2, outer = TRUE, adj = 0)
On this stand, the trees seem to shock concordantly in 2002. But, is this low growth year considered as a shock by the model?
par(mfrow = c(5,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('delta_sck[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=paste0("log(rw) - tree ",t),
display_ylim=c(-12, 1), display_xlim = c(1980, 2010))
abline(v = 2002, lty = 2, col = 'grey40')
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.7, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.4, col="black")
}
mtext(paste0("Stand ", data$uniq_stand_ids[s3]), side = 3, line = -2, outer = TRUE, adj = 0)
Actually, no! The concordance is here captured by the stand-level GP.
par(mfrow = c(1,1), mar = c(2, 5, 3, 1))
names <- sapply(1:31,
function(y) paste0('f_sh[', s3, ',', y, ']'))
util$plot_realizations(samples, names, plot_xs = data$all_years[1:31], N_plots = 50,
main = paste0('Stand ', data$uniq_stand_ids[s]), ylab = 'f_sh',
display_xlim = data$all_years[c(1,31)], display_ylim = c(-2,2))
abline(v = 2002, lty = 2, col = 'grey40')
I think that the shock model is trying only to capture the big shocks (the tail of distribution ?), and thus ignore the medium shocks.
How to improve this?
8 A new implementation to keep only the large shocks
Note
Update Nov. 26 10pm: I try a new method to marginalize only the non-shocks, that does not use the normal-normal conjugancy anymore. This will be useful if we want to use another distribution with a heavier tail for the \delta_\text{shock}.
What if I leave the \delta_\text{shock} explicit in the large shock branch of the mixture model?
I can set: \delta_\text{large shock} \sim \text{normal}(0, {\tau_\text{outer}^\text{shock}}) So, for observations that end up in the non-shocky branch of the model, \delta_\text{large shock} will be inform only by this prior… But then in the generated quantities block I can reconstruct \delta_\text{shock}^\text{post}, which is basically =0 in a non-shock year and \delta_\text{large shock} in a shock year?
Note that in this case I also assume that \tau_\text{inner}^\text{shock} = 0. So I have this in the model block:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalsmallshocks_explicitbigshocks.stan'))[295]) delta_large_sck ~ normal(0, outer_tau_sck);writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalsmallshocks_explicitbigshocks.stan'))[308:318]) for(y in 1:N_years[t]) {
int ys = all_years_idxs_tree[y]-stand_start_years_idxs[s]+1;
log_p0[ys] += log_mix(omega_nonconc_sck,
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]] + f[tree_idxs[y]] + delta_large_sck[tree_idxs[y]], sigma),
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]] + f[tree_idxs[y]], sigma));
log_p1[ys] += log_mix(omega_conc_sck,
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]] + f[tree_idxs[y]] + delta_large_sck[tree_idxs[y]], sigma),
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]] + f[tree_idxs[y]], sigma));
}And in the generated quantities block:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalsmallshocks_explicitbigshocks.stan'))[374:379]) if(sck_state[tree_idxs[y]] == 1) {
delta_sck_post[tree_idxs[y]] = delta_large_sck[tree_idxs[y]];
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]] + delta_sck_post[tree_idxs[y]], sigma);
}else{
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]], sigma);
}I run this model one the small dataset:
data <- readRDS(file = file.path(wd, 'output', 'model', 'data_15nov2025_az592.rds'))
data$N_stand_trees <- array(data$N_stand_trees, dim = 1)
data$N_stand_years <- array(data$N_stand_years, dim = 1)
data$stand_start_years_idxs <- array(data$stand_start_years_idxs, dim = 1)
data$stand_trees_end_idxs <- array(data$stand_trees_end_idxs, dim = 1)
data$stand_trees_start_idxs <- array(data$stand_trees_start_idxs, dim = 1)
if(FALSE){
fit2bis <- stan(file=file.path(wd, 'model', 'stan', 'model9_marginalsmallshocks_explicitbigshocks.stan'),
data=data, seed=5838293, chains = 4,
warmup=1000, iter=2024, refresh=10)
saveRDS(fit2bis, file.path(wd, 'model/shocks/output/marginal', 'fit2bis.rds'))
}else{
fit2bis <- readRDS(file.path(wd, 'model/shocks/output/marginal', 'fit2bis.rds'))
}Diagnostics are not tooootally clean…
diagnostics <- util$extract_hmc_diagnostics(fit2bis)
util$check_all_hmc_diagnostics(diagnostics) Chain 1: E-FMI = 0.197.
Chain 3: E-FMI = 0.169.
Chain 4: E-FMI = 0.185.
E-FMI below 0.2 arise when a funnel-like geometry obstructs how
effectively Hamiltonian trajectories can explore the target
distribution.samples <- util$extract_expectand_vals(fit2bis)
base_samples <- util$filter_expectands(samples,
c('alpha',
'beta_gdd', 'beta_sm', 'beta_vpd',
'rho_sh', 'gamma_sh',
'rho_sp', 'gamma_sp',
'sigma', 'outer_tau_sck',
'phi_sck', 'omega_conc_sck', 'omega_nonconc_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)sigma:
Chain 3: hat{ESS} (42.193) is smaller than desired (100).
outer_tau_sck:
Chain 1: hat{ESS} (58.040) is smaller than desired (100).
Chain 2: hat{ESS} (67.217) is smaller than desired (100).
Chain 3: hat{ESS} (41.556) is smaller than desired (100).
Chain 4: hat{ESS} (37.806) is smaller than desired (100).
omega_nonconc_sck:
Chain 3: hat{ESS} (97.858) is smaller than desired (100).
Small empirical effective sample sizes result in imprecise Markov chain
Monte Carlo estimators.But let’s still give a look at the latent \delta_\text{large shock}:
set.seed(123456)
random_trees <- sample(1:data$N_trees, 4)
par(mfrow = c(2,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('delta_large_sck[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=expression(delta[largeshock]),
display_ylim=c(-10, 10), display_xlim = c(1980, 2010))
}
As expected, their posteriors are informed by the data mostly in 2002 (the main shocky year).
But if we look at the reconstructed \delta_\text{shock}^\text{post}:
par(mfrow = c(2,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('delta_sck_post[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=expression(delta[shock]^post),
display_ylim=c(-12, 1), display_xlim = c(1980, 2010))
}
And the corresponding predictions:
par(mfrow = c(2,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='log(rw)',
display_ylim=c(-12, 1), display_xlim = c(1980, 2010))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.7, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.4, col="black")
}
And the residuals:
par(mfrow = c(2,2), mar = c(2, 5, 3, 1))
for(t in random_trees){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
baseline_values = data$log_rw_obs[idxs_t],
xlab="", ylab=paste0("Residuals - tree ",t),
display_ylim=c(-3, 1), display_xlim = c(1980, 2010),
residual = TRUE)
}
We don’t have the same issue that I had before with the shocks being a bit too large. I think I was doing something wrong with the normal-normal conjugancy? Likely because \sigma is not the variability added to \delta_\text{shock}, but to \mu + f + \delta_\text{shock}.
So when I was doing this:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks_zerotauinner.stan'))[398:406]) if(sck_state[tree_idxs[y]] == 1) {
// we can reconstruct shock posterior using the normal-normal conjugancy
real conjugate_mean = (outer_tau_sck^2 / (outer_tau_sck^2 + sigma^2)) * log_rw_obs[tree_idxs[y]];
real conjugate_sd = sqrt((outer_tau_sck^2 * sigma^2) / (outer_tau_sck^2 + sigma^2));
delta_sck[tree_idxs[y]] = normal_rng(conjugate_mean, conjugate_sd);
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]] + delta_sck[tree_idxs[y]], sigma);
}else{
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]], sigma);
}Maybe the conjugate mean was wrong? I was conditioning on the observation log_rw_obs[tree_idxs[y]], but maybe I should have conditioned on (log_rw_obs[tree_idxs[y]] - mu[y] - f[tree_idxs[y]]), i.e. rather do this:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks_zerotauinner_updated.stan'))[398:407])NA
NA
NA
NA
NA
NA
NA
NA
NA
NAI did not try, though.
9 Exploring the shocks that are not captured
We’ve seen that the shock model is able to capture only the large shocks (missing rings), and let the stand-level GP contort to get the medium ones.
Trying to diagnostic a model that takes hours to run is a pain, so I subset the data to two of the sites I mentioned before: one with large shocks (S14), and one with medium ones (S15).
data <- readRDS(file = file.path(wd, 'output', 'model', 'data_15nov2025_az592_az628.rds'))
if(FALSE){
fit5 <- stan(file=file.path(wd, 'model', 'stan', 'model9_marginalsmallshocks_explicitbigshocks_lessconcordant.stan'),
data=data, seed=5838293, chains = 4,
warmup=1000, iter=2024, refresh=10)
saveRDS(fit5, file.path(wd, 'model/shocks/output/marginal', 'fit5.rds'))
}else{
fit5 <- readRDS(file.path(wd, 'model/shocks/output/marginal', 'fit5.rds'))
}As expected, we also get funnel-related warnings!
diagnostics <- util$extract_hmc_diagnostics(fit5)
util$check_all_hmc_diagnostics(diagnostics) Chain 1: E-FMI = 0.098.
Chain 2: 158 of 1024 transitions (15.43%)
saturated the maximum treedepth of 10.
Chain 2: E-FMI = 0.067.
Chain 3: 134 of 1024 transitions (13.09%)
saturated the maximum treedepth of 10.
Chain 3: E-FMI = 0.096.
Chain 4: 730 of 1024 transitions (71.29%)
saturated the maximum treedepth of 10.
Chain 4: E-FMI = 0.091.
Numerical trajectories that saturate the maximum treedepth have
terminated prematurely. Increasing max_depth above 10 should result in
more expensive, but more efficient, Hamiltonian transitions.
E-FMI below 0.2 arise when a funnel-like geometry obstructs how
effectively Hamiltonian trajectories can explore the target
distribution.samples <- util$extract_expectand_vals(fit5)
base_samples <- util$filter_expectands(samples,
c('alpha',
'beta_gdd', 'beta_sm', 'beta_vpd',
'rho_sh', 'gamma_sh',
'rho_sp', 'gamma_sp',
'sigma', 'outer_tau_sck',
'phi_sck', 'omega_conc_sck', 'omega_nonconc_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)alpha:
Chain 2: hat{ESS} (80.648) is smaller than desired (100).
beta_sm[1]:
Chain 1: hat{ESS} (78.910) is smaller than desired (100).
Chain 3: hat{ESS} (74.330) is smaller than desired (100).
beta_vpd[1]:
Chain 1: hat{ESS} (70.589) is smaller than desired (100).
Chain 3: hat{ESS} (87.303) is smaller than desired (100).
gamma_sh:
Chain 2: hat{ESS} (96.934) is smaller than desired (100).
rho_sp[1]:
Chain 1: hat{ESS} (66.933) is smaller than desired (100).
Chain 2: hat{ESS} (54.624) is smaller than desired (100).
Chain 3: hat{ESS} (77.826) is smaller than desired (100).
sigma:
Split hat{R} (1.133) exceeds 1.1.
Chain 1: hat{ESS} (21.127) is smaller than desired (100).
Chain 2: hat{ESS} (20.674) is smaller than desired (100).
Chain 3: hat{ESS} (22.998) is smaller than desired (100).
Chain 4: hat{ESS} (55.604) is smaller than desired (100).
outer_tau_sck:
Split hat{R} (1.213) exceeds 1.1.
Chain 1: hat{ESS} (13.754) is smaller than desired (100).
Chain 2: hat{ESS} (15.992) is smaller than desired (100).
Chain 3: hat{ESS} (16.284) is smaller than desired (100).
Chain 4: hat{ESS} (19.571) is smaller than desired (100).
phi_sck[2]:
Chain 2: hat{ESS} (42.101) is smaller than desired (100).
Chain 3: hat{ESS} (52.513) is smaller than desired (100).
omega_conc_sck:
Chain 1: hat{ESS} (63.900) is smaller than desired (100).
Chain 2: hat{ESS} (28.831) is smaller than desired (100).
Chain 3: hat{ESS} (58.265) is smaller than desired (100).
omega_nonconc_sck:
Chain 2: hat{ESS} (67.717) is smaller than desired (100).
Split Rhat larger than 1.1 suggests that at least one of the Markov
chains has not reached an equilibrium.
Small empirical effective sample sizes result in imprecise Markov chain
Monte Carlo estimators.Notably, we get \hat{R} warnings for \tau_\text{inner}^\text{shock} and \sigma:
util$plot_pairs_by_chain(samples[['sigma']], expression(sigma), samples[['outer_tau_sck']], expression(tau[shock]^outer))
And quite a lot of issues with our \delta_\text{large shock}:
base_samples <- util$filter_expectands(samples,
c('delta_large_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)delta_large_sck[148]:
Chain 1: Left tail hat{xi} (0.613) exceeds 0.25.
Chain 2: Left tail hat{xi} (0.722) exceeds 0.25.
Chain 3: Left tail hat{xi} (0.603) exceeds 0.25.
Chain 4: Left tail hat{xi} (0.339) exceeds 0.25.
delta_large_sck[149]:
Chain 1: Right tail hat{xi} (0.741) exceeds 0.25.
Chain 1: hat{ESS} (28.717) is smaller than desired (100).
delta_large_sck[151]:
Chain 1: Both left and right tail hat{xi}s (0.667, 0.835) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.532, 0.970) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.466, 0.978) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.796, 1.009) exceed 0.25.
delta_large_sck[152]:
Chain 4: Both left and right tail hat{xi}s (0.463, 0.844) exceed 0.25.
Chain 4: hat{ESS} (26.542) is smaller than desired (100).
delta_large_sck[182]:
Chain 1: Both left and right tail hat{xi}s (0.497, 0.821) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.684, 0.965) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.867, 1.065) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.952, 0.935) exceed 0.25.
delta_large_sck[210]:
Chain 1: Right tail hat{xi} (0.480) exceeds 0.25.
Chain 2: Right tail hat{xi} (0.426) exceeds 0.25.
delta_large_sck[213]:
Chain 3: Left tail hat{xi} (0.416) exceeds 0.25.
delta_large_sck[214]:
Chain 1: Both left and right tail hat{xi}s (0.759, 0.967) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.763, 0.887) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.741, 0.688) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.840, 1.036) exceed 0.25.
delta_large_sck[234]:
Chain 4: Right tail hat{xi} (0.738) exceeds 0.25.
Chain 4: hat{ESS} (24.511) is smaller than desired (100).
delta_large_sck[242]:
Chain 1: Left tail hat{xi} (0.407) exceeds 0.25.
delta_large_sck[244]:
Chain 1: Both left and right tail hat{xi}s (0.744, 0.969) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.639, 0.635) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.852, 0.932) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.781, 0.651) exceed 0.25.
delta_large_sck[276]:
Chain 1: Left tail hat{xi} (0.353) exceeds 0.25.
Chain 2: Left tail hat{xi} (0.532) exceeds 0.25.
Chain 3: Both left and right tail hat{xi}s (0.561, 0.290) exceed 0.25.
Chain 4: Left tail hat{xi} (0.678) exceeds 0.25.
delta_large_sck[334]:
Chain 1: Right tail hat{xi} (0.897) exceeds 0.25.
Chain 2: Both left and right tail hat{xi}s (0.441, 0.577) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.693, 0.642) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.704, 0.930) exceed 0.25.
Chain 3: hat{ESS} (30.083) is smaller than desired (100).
delta_large_sck[337]:
Chain 1: Left tail hat{xi} (0.412) exceeds 0.25.
Chain 2: Both left and right tail hat{xi}s (0.337, 0.452) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.331, 0.736) exceed 0.25.
delta_large_sck[427]:
Chain 1: Left tail hat{xi} (0.810) exceeds 0.25.
Chain 2: Both left and right tail hat{xi}s (0.889, 0.672) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (1.107, 0.672) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (1.110, 0.641) exceed 0.25.
Chain 1: hat{ESS} (70.775) is smaller than desired (100).
delta_large_sck[429]:
Chain 1: Right tail hat{xi} (0.394) exceeds 0.25.
Chain 2: Right tail hat{xi} (0.333) exceeds 0.25.
Chain 3: Right tail hat{xi} (0.459) exceeds 0.25.
Chain 4: Right tail hat{xi} (0.424) exceeds 0.25.
delta_large_sck[430]:
Chain 1: Both left and right tail hat{xi}s (0.906, 0.689) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.662, 0.759) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.943, 0.837) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.687, 0.831) exceed 0.25.
delta_large_sck[520]:
Chain 1: Both left and right tail hat{xi}s (0.803, 0.930) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.726, 1.083) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.779, 0.541) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.925, 0.897) exceed 0.25.
delta_large_sck[543]:
Chain 1: Left tail hat{xi} (0.504) exceeds 0.25.
Chain 2: Left tail hat{xi} (0.583) exceeds 0.25.
Chain 4: Left tail hat{xi} (0.449) exceeds 0.25.
delta_large_sck[551]:
Chain 2: Left tail hat{xi} (0.475) exceeds 0.25.
delta_large_sck[578]:
Chain 1: Left tail hat{xi} (1.035) exceeds 0.25.
Chain 2: Right tail hat{xi} (0.882) exceeds 0.25.
Chain 3: Left tail hat{xi} (1.261) exceeds 0.25.
Chain 4: Both left and right tail hat{xi}s (0.630, 0.382) exceed 0.25.
Chain 1: hat{ESS} (22.459) is smaller than desired (100).
Chain 2: hat{ESS} (22.515) is smaller than desired (100).
Chain 3: hat{ESS} (9.313) is smaller than desired (100).
delta_large_sck[585]:
Chain 1: Left tail hat{xi} (0.262) exceeds 0.25.
Chain 2: Both left and right tail hat{xi}s (0.512, 0.327) exceed 0.25.
Chain 3: Left tail hat{xi} (0.376) exceeds 0.25.
Chain 4: Left tail hat{xi} (0.514) exceeds 0.25.
delta_large_sck[604]:
Chain 1: Both left and right tail hat{xi}s (0.849, 0.742) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.613, 0.714) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.765, 0.636) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.666, 0.596) exceed 0.25.
delta_large_sck[650]:
Chain 1: Both left and right tail hat{xi}s (0.513, 0.921) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.552, 0.881) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.840, 1.059) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.769, 1.043) exceed 0.25.
delta_large_sck[667]:
Chain 3: Right tail hat{xi} (0.298) exceeds 0.25.
Chain 4: Left tail hat{xi} (0.701) exceeds 0.25.
Chain 4: hat{ESS} (61.346) is smaller than desired (100).
delta_large_sck[675]:
Chain 1: Both left and right tail hat{xi}s (0.712, 0.967) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.837, 0.983) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.736, 1.050) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.834, 0.399) exceed 0.25.
delta_large_sck[679]:
Chain 1: Left tail hat{xi} (1.535) exceeds 0.25.
Chain 3: Both left and right tail hat{xi}s (0.580, 0.513) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.352, 0.763) exceed 0.25.
Split hat{R} (1.159) exceeds 1.1.
Chain 1: hat{ESS} (7.323) is smaller than desired (100).
Chain 4: hat{ESS} (84.180) is smaller than desired (100).
delta_large_sck[680]:
Chain 1: Left tail hat{xi} (0.352) exceeds 0.25.
Chain 3: Left tail hat{xi} (0.256) exceeds 0.25.
delta_large_sck[687]:
Chain 1: Right tail hat{xi} (0.273) exceeds 0.25.
delta_large_sck[691]:
Chain 1: Both left and right tail hat{xi}s (0.910, 0.635) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (1.070, 0.635) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.698, 0.799) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.418, 0.580) exceed 0.25.
delta_large_sck[752]:
Chain 2: Right tail hat{xi} (0.507) exceeds 0.25.
Chain 3: Both left and right tail hat{xi}s (0.326, 1.271) exceed 0.25.
Chain 4: Right tail hat{xi} (1.095) exceeds 0.25.
Chain 2: hat{ESS} (55.293) is smaller than desired (100).
Chain 3: hat{ESS} (15.190) is smaller than desired (100).
Chain 4: hat{ESS} (31.414) is smaller than desired (100).
delta_large_sck[848]:
Chain 1: Left tail hat{xi} (0.511) exceeds 0.25.
Chain 2: Left tail hat{xi} (0.417) exceeds 0.25.
delta_large_sck[967]:
Chain 1: Both left and right tail hat{xi}s (0.540, 0.720) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.616, 0.980) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.739, 0.989) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.796, 0.845) exceed 0.25.
delta_large_sck[977]:
Chain 1: Both left and right tail hat{xi}s (0.666, 0.967) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.780, 1.017) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.749, 0.283) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.751, 0.418) exceed 0.25.
delta_large_sck[1020]:
Chain 1: Both left and right tail hat{xi}s (0.827, 0.880) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.754, 0.716) exceed 0.25.
Chain 3: Left tail hat{xi} (0.401) exceeds 0.25.
Chain 4: Both left and right tail hat{xi}s (0.510, 0.285) exceed 0.25.
delta_large_sck[1063]:
Chain 1: Right tail hat{xi} (0.430) exceeds 0.25.
Chain 2: Right tail hat{xi} (0.454) exceeds 0.25.
delta_large_sck[1118]:
Chain 1: Both left and right tail hat{xi}s (0.348, 0.599) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.978, 0.856) exceed 0.25.
delta_large_sck[1191]:
Chain 1: Both left and right tail hat{xi}s (0.754, 0.485) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.729, 0.596) exceed 0.25.
Chain 3: Left tail hat{xi} (0.832) exceeds 0.25.
Chain 4: Left tail hat{xi} (0.764) exceeds 0.25.
delta_large_sck[1200]:
Chain 3: Left tail hat{xi} (0.888) exceeds 0.25.
Chain 3: hat{ESS} (20.927) is smaller than desired (100).
delta_large_sck[1317]:
Chain 1: Both left and right tail hat{xi}s (0.787, 0.356) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.559, 0.289) exceed 0.25.
Chain 3: Right tail hat{xi} (0.343) exceeds 0.25.
Chain 4: Both left and right tail hat{xi}s (0.647, 0.669) exceed 0.25.
Chain 4: hat{ESS} (60.554) is smaller than desired (100).
delta_large_sck[1321]:
Chain 1: Left tail hat{xi} (0.660) exceeds 0.25.
Chain 2: Both left and right tail hat{xi}s (0.993, 0.540) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (1.217, 0.755) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.888, 0.879) exceed 0.25.
delta_large_sck[1325]:
Chain 1: Left tail hat{xi} (0.277) exceeds 0.25.
Chain 2: Left tail hat{xi} (0.402) exceeds 0.25.
Chain 3: Left tail hat{xi} (0.461) exceeds 0.25.
delta_large_sck[1450]:
Chain 1: Left tail hat{xi} (0.530) exceeds 0.25.
Chain 2: Left tail hat{xi} (0.307) exceeds 0.25.
delta_large_sck[1462]:
Chain 2: Left tail hat{xi} (0.254) exceeds 0.25.
delta_large_sck[1483]:
Chain 1: Right tail hat{xi} (0.457) exceeds 0.25.
Chain 2: Right tail hat{xi} (0.545) exceeds 0.25.
delta_large_sck[1505]:
Chain 1: Both left and right tail hat{xi}s (0.500, 0.890) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.884, 0.726) exceed 0.25.
Chain 3: Right tail hat{xi} (0.304) exceeds 0.25.
delta_large_sck[1536]:
Chain 1: Right tail hat{xi} (0.272) exceeds 0.25.
delta_large_sck[1579]:
Chain 1: Both left and right tail hat{xi}s (0.746, 0.891) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.842, 0.981) exceed 0.25.
Chain 3: Left tail hat{xi} (0.532) exceeds 0.25.
Chain 4: Both left and right tail hat{xi}s (0.827, 1.023) exceed 0.25.
Chain 3: hat{ESS} (35.825) is smaller than desired (100).
delta_large_sck[1665]:
Chain 1: Both left and right tail hat{xi}s (0.699, 0.429) exceed 0.25.
Chain 2: Left tail hat{xi} (0.720) exceeds 0.25.
Chain 4: Left tail hat{xi} (0.494) exceeds 0.25.
delta_large_sck[1708]:
Chain 1: Both left and right tail hat{xi}s (0.868, 0.621) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.994, 0.600) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (1.179, 0.691) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (1.133, 0.842) exceed 0.25.
Chain 3: hat{ESS} (79.090) is smaller than desired (100).
delta_large_sck[1712]:
Chain 1: Both left and right tail hat{xi}s (0.576, 0.585) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.344, 0.584) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.389, 0.736) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.705, 0.813) exceed 0.25.
delta_large_sck[1713]:
Chain 1: Left tail hat{xi} (0.412) exceeds 0.25.
Chain 2: Left tail hat{xi} (0.453) exceeds 0.25.
Chain 3: Left tail hat{xi} (0.596) exceeds 0.25.
Chain 4: Left tail hat{xi} (0.540) exceeds 0.25.
delta_large_sck[1741]:
Chain 1: Both left and right tail hat{xi}s (0.680, 0.881) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (1.165, 0.914) exceed 0.25.
Chain 3: Right tail hat{xi} (0.801) exceeds 0.25.
Chain 4: Right tail hat{xi} (0.498) exceeds 0.25.
delta_large_sck[1795]:
Chain 1: Both left and right tail hat{xi}s (0.700, 1.163) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.572, 1.100) exceed 0.25.
Chain 3: Both left and right tail hat{xi}s (0.541, 0.835) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.715, 1.092) exceed 0.25.
Chain 1: hat{ESS} (35.306) is smaller than desired (100).
delta_large_sck[1880]:
Chain 1: Both left and right tail hat{xi}s (0.310, 1.435) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.720, 1.040) exceed 0.25.
Chain 3: Left tail hat{xi} (0.626) exceeds 0.25.
Chain 4: Both left and right tail hat{xi}s (0.821, 1.135) exceed 0.25.
Chain 1: hat{ESS} (18.884) is smaller than desired (100).
Chain 3: hat{ESS} (14.217) is smaller than desired (100).
delta_large_sck[1892]:
Chain 2: Both left and right tail hat{xi}s (0.333, 0.685) exceed 0.25.
delta_large_sck[1913]:
Chain 1: Both left and right tail hat{xi}s (0.811, 0.917) exceed 0.25.
Chain 2: Both left and right tail hat{xi}s (0.864, 0.843) exceed 0.25.
delta_large_sck[1923]:
Chain 1: Left tail hat{xi} (0.254) exceeds 0.25.
Chain 2: Both left and right tail hat{xi}s (0.480, 0.870) exceed 0.25.
Chain 4: Both left and right tail hat{xi}s (0.581, 0.818) exceed 0.25.
Chain 1: hat{ESS} (71.636) is smaller than desired (100).
Chain 3: hat{ESS} (11.372) is smaller than desired (100).
Chain 4: hat{ESS} (73.218) is smaller than desired (100).
Large tail hat{xi}s suggest that the expectand might not be
sufficiently integrable.
Split Rhat larger than 1.1 suggests that at least one of the Markov
chains has not reached an equilibrium.
Small empirical effective sample sizes result in imprecise Markov chain
Monte Carlo estimators.Let’s look at one of the shock that doesn’t have a nice geometry:
util$plot_div_pairs(paste0('delta_large_sck[', 1712, ']'), 'sigma', samples, diagnostics, transforms = list('sigma' = 1))
This shock corresponds to the dashed line below:
par(mfrow = c(3,1))
t <- 46
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='log_rw_pred',
display_ylim=c(-2, 2), display_xlim = c(1980, 2021))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=1, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.6, col="black")
abline(v = data$years[1712], lty = 2)
names <- sapply(idxs_t,
function(x) paste0('delta_sck_post[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='delta_sck_post',
display_ylim=c(-2, 2), display_xlim = c(1980, 2021))
abline(v = data$years[1712], lty = 2)
names <- sapply(1:43,
function(x) paste0('f_sh[', data$stand_idxs[t],',',x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='Stand GP',
display_ylim=c(-2, 2), display_xlim = c(1980, 2021))
abline(v = data$years[1712], lty = 2)
As a side note, the diagnostics are cleaner when we keep everything marginalized and we reconstruct the \delta_\text{shock} with the normal-normal conjugancy (corrected):
data <- readRDS(file = file.path(wd, 'output', 'model', 'data_15nov2025_az592_az628.rds'))
if(FALSE){
fit6 <- stan(file=file.path(wd, 'model', 'stan', 'model9_marginalshocks_zerotauinner_updated.stan'),
data=data, seed=5838293, chains = 4,
warmup=1000, iter=2024, refresh=10)
saveRDS(fit6, file.path(wd, 'model/shocks/output/marginal', 'fit6.rds'))
}else{
fit6 <- readRDS(file.path(wd, 'model/shocks/output/marginal', 'fit6.rds'))
}diagnostics <- util$extract_hmc_diagnostics(fit6)
util$check_all_hmc_diagnostics(diagnostics) All Hamiltonian Monte Carlo diagnostics are consistent with reliable
Markov chain Monte Carlo.samples <- util$extract_expectand_vals(fit6)
base_samples <- util$filter_expectands(samples,
c('alpha',
'beta_gdd', 'beta_sm', 'beta_vpd',
'rho_sh', 'gamma_sh',
'rho_sp', 'gamma_sp',
'sigma', 'outer_tau_sck',
'phi_sck', 'omega_conc_sck', 'omega_nonconc_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)beta_sm[1]:
Chain 4: hat{ESS} (68.836) is smaller than desired (100).
beta_vpd[1]:
Chain 4: hat{ESS} (73.587) is smaller than desired (100).
rho_sp[1]:
Chain 2: hat{ESS} (57.838) is smaller than desired (100).
Small empirical effective sample sizes result in imprecise Markov chain
Monte Carlo estimators.And for the tree mentioned above we obtain:
par(mfrow = c(3,1))
t <- 46
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='log_rw_pred',
display_ylim=c(-2, 2), display_xlim = c(1980, 2021))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=1, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.6, col="black")
abline(v = data$years[1712], lty = 2)
names <- sapply(idxs_t,
function(x) paste0('delta_sck[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='delta_sck',
display_ylim=c(-2, 2), display_xlim = c(1980, 2021))
abline(v = data$years[1712], lty = 2)
names <- sapply(1:43,
function(x) paste0('f_sh[', data$stand_idxs[t],',',x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='Stand GP',
display_ylim=c(-2, 2), display_xlim = c(1980, 2021))
abline(v = data$years[1712], lty = 2)
10 Move to a Cauchy?
Note
Update Dec. 1 11am: Reading a bit more about heavy-tailed distribution, I wonder if I should rather consider a Student-t distribution? Because the marginalization seems kind of easier? And with nu = 1, we obtain a Cauchy, right? I will likely try in a latter section…
We might need a heavier-tail distribution to catch these heterogeneous shocks!
10.1 Pseudo-code
Before any marginalization, we had this:
target += normal_lpdf(y[n] | mu[n] + delta_sck[n], sigma)
target += normal_lpdf(delta_sck[n] | 0, outer_tau_sck)Moving to a Cauchy, we would have this:
target += normal_lpdf(y[n] | mu[n] + delta_sck[n], sigma)
target += cauchy_lpdf(delta_sck[n] | 0, outer_tau_sck)\tau_\text{outer}^\text{shock} is now the scale of the Cauchy distribution: \delta_\text{shock} \sim \text{Cauchy}(0, \tau_\text{outer}^\text{shock}) I came across Mike’s writing on Cauchy, and in particular the second implementation. If I understood correctly, the Cauchy is a scale mixture of Gaussian density functions, and in particular we would have: \delta_\text{shock} \mid \tau_\text{outer}^\text{aux} \sim \text{normal}(0, \tau_\text{outer}^\text{aux}) {\tau_\text{outer}^\text{aux}}^2 \sim \text{InvGamma}(0.5, 0.5*{\tau_\text{outer}^\text{shock}}^2) So, I think we could write:
target += normal_lpdf(y[n] | mu[n] + delta_sck[n], sigma)
target += normal_lpdf(delta_sck[n] | 0, sqrt(outer_tau2_aux))
target += inv_gamma_lpdf(outer_tau2_aux | 0.5, 0.5*outer_tau_sck^2)And thus marginalize the two first lines and obtain:
target += normal_lpdf(y[n] | mu[n], sqrt(outer_tau2_aux+sigma^2))
target += inv_gamma_lpdf(outer_tau2_aux | 0.5, 0.5*outer_tau_sck^2)Which may allow to keep the normal-normal conjugancy to reconstruct \delta_\text{shock} afterwards..? I never thought I would be so much into conjugate distributions.
In terms of domain expertise, we have to rethink about our prior on \tau_\text{outer}^\text{shock}. Playing with this magic tool, we could set a prior to have 0\lesssim\tau_\text{outer}^\text{shock}\lesssim0.4, as it would imply roughly shocks between -3 and 3 (Q5-Q95), and large shocks between -12 and 12 (Q1-Q99). So something like:
outer_tau_sck ~ normal(0, 0.4/ 2.57);10.2 Actual Stan code
Our observational model is now:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks_zerotauinner_cauchy.stan'))[281:318]) // Cauchy is a scale mixture of Gaussian density functions!
outer_tau_sck ~ normal(0, 0.2/ 2.57); // outer_tau_sck is now the scale of Cauchy. With outer_tau_sck = 1, we get a Cauchy between -12 and 12
outer_tau2_aux ~ inv_gamma(0.5, 0.5 * square(outer_tau_sck));
// Observational model with marginalized shocks
for (s in 1:N_stands) {
array[N_stand_trees[s]] int stand_trees_idxs = linspaced_int_array(N_stand_trees[s],
stand_trees_start_idxs[s], stand_trees_end_idxs[s]);
vector[N_stand_years[s]] log_p0 = rep_vector(0, N_stand_years[s]);
vector[N_stand_years[s]] log_p1 = rep_vector(0, N_stand_years[s]);
for(ts in 1:N_stand_trees[s]){
int t = stand_trees_idxs[ts];
array[N_years[t]] int tree_idxs = linspaced_int_array(N_years[t], tree_start_idxs[t], tree_end_idxs[t]);
array[N_years[t]] int all_years_idxs_tree = all_years_idxs[tree_idxs];
for(y in 1:N_years[t]) {
int ys = all_years_idxs_tree[y]-stand_start_years_idxs[s]+1;
log_p0[ys] += log_mix(omega_nonconc_sck,
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sqrt(outer_tau2_aux[tree_idxs[y]] + sigma^2)),
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sigma));
log_p1[ys] += log_mix(omega_conc_sck,
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sqrt(outer_tau2_aux[tree_idxs[y]] + sigma^2)),
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sigma));
}
}
for(y in 1:N_stand_years[s]) {
target += log_mix(phi_sck[s], log_p1[y], log_p0[y]);
}
}With the following generated quantities:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks_zerotauinner_cauchy.stan'))[357:416])generated quantities {
vector[N] delta_sck = rep_vector(0,N); // latent amplitude of shock
array[N] int sck_state; // latent state, zt = 0 or zt = 1
array[N] real log_rw_pred;
for (t in 1:N_trees) {
array[N_years[t]] int tree_idxs
= linspaced_int_array(N_years[t],
tree_start_idxs[t],
tree_end_idxs[t]);
array[N_years[t]] int all_years_idxs_tree = all_years_idxs[tree_idxs];
array[N_years[t]] real years_tree = years[tree_idxs];
int stand_idx = stand_idxs[t];
int species_idx = species_idxs[t];
vector[N_years[t]] gdd_obs_tree = gdd_obs[tree_idxs];
vector[N_years[t]] sm_obs_tree = sm_obs[tree_idxs];
vector[N_years[t]] vpd_obs_tree = vpd_obs[tree_idxs];
vector[N_years[t]] mu = alpha
+ beta_gdd[species_idx] * (gdd_obs_tree - gdd0)
+ beta_sm[species_idx] * (sm_obs_tree - sm0)
+ beta_vpd[species_idx] * (vpd_obs_tree - vpd0)
+ f_sh[stand_idx, all_years_idxs_tree];
// mixture weight for shock
real mw_shock = phi_sck[stand_idx]*omega_conc_sck + (1-phi_sck[stand_idx])*omega_nonconc_sck;
for(y in 1:N_years[t]){
// mixture components
real log_pshock = log(mw_shock) + normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[y] + f[tree_idxs[y]], sqrt(outer_tau2_aux[tree_idxs[y]] + sigma^2));
real log_pshock_plus_pnonshock = log_mix(mw_shock,
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[y] + f[tree_idxs[y]], sqrt(outer_tau2_aux[tree_idxs[y]] + sigma^2)),
normal_lpdf(log_rw_obs[tree_idxs[y]] | mu[y] + f[tree_idxs[y]], sigma));
// probability of shock state
real lambda_shock = exp(log_pshock - log_pshock_plus_pnonshock);
sck_state[tree_idxs[y]] = bernoulli_rng(lambda_shock); // or something like categorical_rng(lambda_shock);?
if(sck_state[tree_idxs[y]] == 1) {
// we can reconstruct shock posterior using the normal-normal conjugancy
real residual = log_rw_obs[tree_idxs[y]] - mu[y] - f[tree_idxs[y]];
real conjugate_mean = (outer_tau2_aux[tree_idxs[y]] / (outer_tau2_aux[tree_idxs[y]] + sigma^2)) * residual;
real conjugate_sd = sqrt((outer_tau2_aux[tree_idxs[y]] * sigma^2) / (outer_tau2_aux[tree_idxs[y]] + sigma^2));
delta_sck[tree_idxs[y]] = normal_rng(conjugate_mean, conjugate_sd);
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]] + delta_sck[tree_idxs[y]], sigma);
}else{
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]], sigma);
}
}
}
}So, let’s run this!
No, wait…
11 Time to deal with those zeros
Before digging too much on what can make the model not to capture intermediate and large shocks together, I should first consider to better define what is a large shock. In our two sites example, most large shocks actually correspond to missing rings—or at least below the measurement precision, which should be around 1\mathrm{e}{-3} ~\text{mm}! And until now we were considering something like: \log(y + \epsilon) ~ \text{with} ~ \epsilon = 1\mathrm{e}{-4} ~\text{mm} Boooo! Loosers.
We should rather consider that the probability of having a missing ring corresponds to the probability of getting values from -Inf to some threshold epsilon (which can be fixed, or a parameter to infer).
if(rw_obs[tree_idxs[y]] > epsilon){
target += normal_lpdf(log(rw_obs[tree_idxs[y]]) | mu[n], sqrt(outer_tau2_aux+sigma^2))
}else{
target += normal_lcdf(log(epsilon)) | mu[n], sqrt(outer_tau2_aux+sigma^2))
}11.1 A model with censoring
We simply modify the model to have this:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks_zerotauinner_cauchy_censoring.stan'))[315:339])
if(rw_obs[tree_idxs[y]] > epsilon){
log_p0[ys] += log_mix(omega_nonconc_sck,
normal_lpdf(log(rw_obs[tree_idxs[y]]) | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sqrt(outer_tau2_aux[tree_idxs[y]] + sigma^2)),
normal_lpdf(log(rw_obs[tree_idxs[y]]) | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sigma));
log_p1[ys] += log_mix(omega_conc_sck,
normal_lpdf(log(rw_obs[tree_idxs[y]]) | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sqrt(outer_tau2_aux[tree_idxs[y]] + sigma^2)),
normal_lpdf(log(rw_obs[tree_idxs[y]]) | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sigma));
}else{
log_p0[ys] += log_mix(omega_nonconc_sck,
normal_lcdf(log(epsilon)| mu[tree_idxs[y]]
+ f[tree_idxs[y]], sqrt(outer_tau2_aux[tree_idxs[y]] + sigma^2)),
normal_lcdf(log(epsilon) | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sigma));
log_p1[ys] += log_mix(omega_conc_sck,
normal_lcdf(log(epsilon) | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sqrt(outer_tau2_aux[tree_idxs[y]] + sigma^2)),
normal_lcdf(log(epsilon) | mu[tree_idxs[y]]
+ f[tree_idxs[y]], sigma));
}But I am less sure on how to handle this in the generated quantities block.
For now I did:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks_zerotauinner_cauchy_censoring.stan'))[442:466])
if(rw_obs[tree_idxs[y]] > epsilon){
if(sck_state[tree_idxs[y]] == 1) {
// we can reconstruct shock posterior using the normal-normal conjugancy
real residual = log(rw_obs[tree_idxs[y]]) - mu[y] - f[tree_idxs[y]];
real conjugate_mean = (outer_tau2_aux[tree_idxs[y]] / (outer_tau2_aux[tree_idxs[y]] + sigma^2)) * residual;
real conjugate_sd = sqrt((outer_tau2_aux[tree_idxs[y]] * sigma^2) / (outer_tau2_aux[tree_idxs[y]] + sigma^2));
delta_sck[tree_idxs[y]] = normal_rng(conjugate_mean, conjugate_sd);
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]] + delta_sck[tree_idxs[y]], sigma);
}else{
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]], sigma);
}
}else{
if(sck_state[tree_idxs[y]] == 1) {
// sample from a truncated normal distribution? between -inf and log(epsilon)
real log_y = normal_ub_rng(mu[y] + f[tree_idxs[y]], sqrt(outer_tau2_aux[tree_idxs[y]] + sigma^2), log(epsilon));
real residual = log_y - mu[y] - f[tree_idxs[y]];
real conjugate_mean = (outer_tau2_aux[tree_idxs[y]] / (outer_tau2_aux[tree_idxs[y]] + sigma^2)) * residual;
real conjugate_sd = sqrt((outer_tau2_aux[tree_idxs[y]] * sigma^2) / (outer_tau2_aux[tree_idxs[y]] + sigma^2));
delta_sck[tree_idxs[y]] = normal_rng(conjugate_mean, conjugate_sd);
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]] + delta_sck[tree_idxs[y]], sigma);
}else{
log_rw_pred[tree_idxs[y]] = normal_rng(mu[y] + f[tree_idxs[y]], sigma);
}
}And based on Stan manual I defined:
writeLines(
readLines(file.path(wd, 'model', 'stan',
'model9_marginalshocks_zerotauinner_cauchy_censoring.stan'))[57:67]) // see truncated-random-number-generation on Stan's user guide
real normal_ub_rng(real mu, real sigma, real ub) {
real p_ub = normal_cdf(ub, mu, sigma);
// deal with the limit case of the truncated normal
// to avoid error: "Exception: uniform_rng: Upper bound parameter is 0"
if (p_ub == 0)
return ub;
real u = uniform_rng(0, p_ub);
real y = mu + sigma * inv_Phi(u);
return y;
}Currently I just set: real epsilon = 1e-3;.
11.2 Test on one shocky site
data <- readRDS(file = file.path(wd, 'output', 'model', 'data_15nov2025_az592.rds'))
data$N_stand_trees <- array(data$N_stand_trees, dim = 1)
data$N_stand_years <- array(data$N_stand_years, dim = 1)
data$stand_start_years_idxs <- array(data$stand_start_years_idxs, dim = 1)
data$stand_trees_end_idxs <- array(data$stand_trees_end_idxs, dim = 1)
data$stand_trees_start_idxs <- array(data$stand_trees_start_idxs, dim = 1)
data$rw_obs <- ifelse(data$log_rw_obs == log(1e-4), 0, exp(data$log_rw_obs))
if(FALSE){
fit <- stan(file=file.path(wd, 'model', 'stan', 'model9_marginalshocks_zerotauinner_cauchy_censoring.stan'),
data=data, seed=5838293, chains = 4,
warmup=1000, iter=2024, refresh=10)
saveRDS(fit, file.path(wd, 'model/shocks/output/marginal', 'fit_cauchy1.rds'))
}else{
fit <- readRDS(file.path(wd, 'model/shocks/output/marginal', 'fit_cauchy1.rds'))
}Diagnostics are fine!
diagnostics <- util$extract_hmc_diagnostics(fit)
util$check_all_hmc_diagnostics(diagnostics) All Hamiltonian Monte Carlo diagnostics are consistent with reliable
Markov chain Monte Carlo.samples <- util$extract_expectand_vals(fit)
base_samples <- util$filter_expectands(samples,
c('alpha',
'beta_gdd', 'beta_sm', 'beta_vpd',
'rho_sh', 'gamma_sh',
'rho_sp', 'gamma_sp',
'sigma', 'outer_tau_sck',
'phi_sck', 'omega_conc_sck', 'omega_nonconc_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)outer_tau_sck:
Chain 1: hat{ESS} (70.634) is smaller than desired (100).
Chain 2: hat{ESS} (61.068) is smaller than desired (100).
Small empirical effective sample sizes result in imprecise Markov chain
Monte Carlo estimators.The shock model now seems to capture the different amplitudes of the 2002 shock across trees. Notice the nice prediction uncertainty on the large shocks.
I wonder if our domain expertise should make us consider 1996 as a shock? Right now it is captured by the stand-level GP, with a little additional shock for WMP03.
par(mfcol = c(4,2))
for(t in 1:4){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='log_rw_pred', main = paste0('Tree ', data$uniq_tree_ids[t]),
display_ylim=c(-10, 2), display_xlim = c(1980, 2010))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=1, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.6, col="black")
abline(v = c(2002,1996), lty = 3, col="darkgrey")
names <- sapply(idxs_t,
function(x) paste0('sck_state[', x, ']'))
ticks <- c(unlist(lapply(seq(1980, 2000,10), function(x) c(x, rep(NA, 9)))), '2010')
util$plot_disc_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='Shock state',
display_ylim=c(-0.1, 1.1),
xticklabs = ticks)
abline(v = c(2002,1996), lty = 3, col="darkgrey")
names <- sapply(idxs_t,
function(x) paste0('delta_sck[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=expression(delta[shock]),
display_ylim=c(-10, 2), display_xlim = c(1980, 2010))
abline(v = c(2002,1996), lty = 3, col="darkgrey")
names <- sapply(1:data$N_all_years,
function(x) paste0('f_sh[', data$stand_idxs[t],',',x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$all_years,
xlab="", ylab='Stand GP',
display_ylim=c(-2, 2), display_xlim = c(1980, 2010))
abline(v = c(2002,1996), lty = 3, col="darkgrey")
}
And the parameters related to the shocks:
# Quickly look at posteriors
par(cex.lab=1.1)
nf <- layout(
matrix(c(1,1,1,2,2,2,
3,3,4,4,5,5),
ncol=6, byrow=TRUE)
)
util$plot_expectand_pushforward(samples[['sigma']], 50,
flim = c(0,0.5),
display_name=expression(sigma))
xs <- seq(0, 2, 0.01)
ys <- dnorm(xs, 0, log(1.1) / 2.57)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples[['outer_tau_sck']], 50,
flim = c(0,1),
display_name=expression(tau[shock]^outer))
xs <- seq(0, 1, 0.01)
ys <- dnorm(xs, 0, 0.4 / 2.57)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples[['phi_sck[1]']], 50,
flim = c(0,1),
display_name=expression(phi[schock]))
xs <- seq(0, 1, 0.01)
ys <- dbeta(xs, 2, 20)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples[['omega_conc_sck']], 50,
flim = c(0.5,1),
display_name=expression(omega[shock]^concordant))
xs <- seq(0.5, 1, 0.01)
ys <- dbeta(xs, 230, 14)
lines(xs, ys, lwd=2.5, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples[['omega_nonconc_sck']], 50,
flim = c(0,0.5),
display_name=expression(omega[shock]^nonconcordant))
xs <- seq(0, 1, 0.01)
ys <- dbeta(xs, 1, 20)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
11.3 Test on two sites
11.3.1 Posterior quantification
data <- readRDS(file = file.path(wd, 'output', 'model', 'data_15nov2025_az592_az628.rds'))
data$rw_obs <- ifelse(data$log_rw_obs == log(1e-4), 0, exp(data$log_rw_obs))
if(FALSE){
fit <- stan(file=file.path(wd, 'model', 'stan', 'model9_marginalshocks_zerotauinner_cauchy_censoring.stan'),
data=data, seed=5838293, chains = 4,
warmup=1000, iter=2024, refresh=10)
saveRDS(fit, file.path(wd, 'model/shocks/output/marginal', 'fit_cauchy2.rds'))
}else{
fit <- readRDS(file.path(wd, 'model/shocks/output/marginal', 'fit_cauchy2.rds'))
}Diagnostics are not very clean…
diagnostics <- util$extract_hmc_diagnostics(fit)
util$check_all_hmc_diagnostics(diagnostics) Chain 1: 596 of 1024 transitions (58.20%)
saturated the maximum treedepth of 10.
Chain 3: 653 of 1024 transitions (63.77%)
saturated the maximum treedepth of 10.
Numerical trajectories that saturate the maximum treedepth have
terminated prematurely. Increasing max_depth above 10 should result in
more expensive, but more efficient, Hamiltonian transitions.samples2 <- util$extract_expectand_vals(fit)
base_samples <- util$filter_expectands(samples2,
c('alpha',
'beta_gdd', 'beta_sm', 'beta_vpd',
'rho_sh', 'gamma_sh',
'rho_sp', 'gamma_sp',
'sigma', 'outer_tau_sck',
'phi_sck', 'omega_conc_sck', 'omega_nonconc_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)alpha:
Chain 2: hat{ESS} (83.172) is smaller than desired (100).
beta_gdd[1]:
Chain 2: hat{ESS} (71.944) is smaller than desired (100).
Chain 4: hat{ESS} (80.468) is smaller than desired (100).
beta_sm[1]:
Chain 2: hat{ESS} (75.391) is smaller than desired (100).
Chain 4: hat{ESS} (83.391) is smaller than desired (100).
beta_vpd[1]:
Chain 2: hat{ESS} (74.092) is smaller than desired (100).
Chain 4: hat{ESS} (79.809) is smaller than desired (100).
gamma_sh:
Chain 4: hat{ESS} (94.180) is smaller than desired (100).
rho_sp[1]:
Chain 4: hat{ESS} (91.387) is smaller than desired (100).
outer_tau_sck:
Chain 2: hat{ESS} (59.827) is smaller than desired (100).
Chain 4: hat{ESS} (57.458) is smaller than desired (100).
Small empirical effective sample sizes result in imprecise Markov chain
Monte Carlo estimators.11.3.2 Stand 1
Unfortunately, 2002 is not considered as a shock for the tree WMP01 anymore…
par(mfcol = c(4,2))
for(t in 1:4){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples2, names, data$years[idxs_t],
xlab="", ylab='log_rw_pred', main = paste0('Tree ', data$uniq_tree_ids[t]),
display_ylim=c(-10, 2), display_xlim = c(1980, 2010))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=1, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.6, col="black")
abline(v = c(2002,1996), lty = 3, col="darkgrey")
names <- sapply(idxs_t,
function(x) paste0('sck_state[', x, ']'))
ticks <- c(unlist(lapply(seq(1980, 2000,10), function(x) c(x, rep(NA, 9)))), '2010')
util$plot_disc_pushforward_quantiles(samples2, names, data$years[idxs_t],
xlab="", ylab='Shock state',
display_ylim=c(-0.1, 1.1),
xticklabs = ticks)
abline(v = c(2002,1996), lty = 3, col="darkgrey")
names <- sapply(idxs_t,
function(x) paste0('delta_sck[', x, ']'))
util$plot_conn_pushforward_quantiles(samples2, names, data$years[idxs_t],
xlab="", ylab=expression(delta[shock]),
display_ylim=c(-10, 2), display_xlim = c(1980, 2010))
abline(v = c(2002,1996), lty = 3, col="darkgrey")
names <- sapply(1:data$N_all_years,
function(x) paste0('f_sh[', data$stand_idxs[t],',',x, ']'))
util$plot_conn_pushforward_quantiles(samples2, names, data$all_years,
xlab="", ylab='Stand GP',
display_ylim=c(-2, 2), display_xlim = c(1980, 2010))
abline(v = c(2002,1996), lty = 3, col="darkgrey")
}
11.3.3 Stand 2
In this stand, we would also expect 2002 to be considered as a shock… and it’s not.
par(mfcol = c(4,2))
for(t in 41:44){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples2, names, data$years[idxs_t],
xlab="", ylab='log_rw_pred', main = paste0('Tree ', data$uniq_tree_ids[t]),
display_ylim=c(-10, 2), display_xlim = c(1980, 2022))
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=1, col="white")
points(data$years[idxs_t], data$log_rw_obs[idxs_t], pch=16, cex=0.6, col="black")
abline(v = c(2002,2006), lty = 3, col="darkgrey")
names <- sapply(idxs_t,
function(x) paste0('sck_state[', x, ']'))
ticks <- c(unlist(lapply(seq(1980, 2010,10), function(x) c(x, rep(NA, 9)))), '2020',NA,NA)
util$plot_disc_pushforward_quantiles(samples2, names, data$years[idxs_t],
xlab="", ylab='Shock state',
display_ylim=c(-0.1, 1.1),
xticklabs = ticks)
abline(v = c(2002,2006), lty = 3, col="darkgrey")
names <- sapply(idxs_t,
function(x) paste0('delta_sck[', x, ']'))
util$plot_conn_pushforward_quantiles(samples2, names, data$years[idxs_t],
xlab="", ylab=expression(delta[shock]),
display_ylim=c(-10, 2), display_xlim = c(1980, 2022))
abline(v = c(2002,2006), lty = 3, col="darkgrey")
names <- sapply(1:data$N_all_years,
function(x) paste0('f_sh[', data$stand_idxs[t],',',x, ']'))
util$plot_conn_pushforward_quantiles(samples2, names, data$all_years,
xlab="", ylab='Stand GP',
display_ylim=c(-2, 2), display_xlim = c(1980, 2022))
abline(v = c(2002,2006), lty = 3, col="darkgrey")
}
11.3.4 Posterior inferences
# Quickly look at posteriors
par(cex.lab=1.1)
nf <- layout(
matrix(c(1,1,1,2,2,2,
3,3,4,4,5,5),
ncol=6, byrow=TRUE)
)
util$plot_expectand_pushforward(samples2[['sigma']], 50,
flim = c(0,0.5),
display_name=expression(sigma))
xs <- seq(0, 2, 0.01)
ys <- dnorm(xs, 0, log(1.1) / 2.57)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples2[['outer_tau_sck']], 50,
flim = c(0,1),
display_name=expression(tau[shock]^outer))
xs <- seq(0, 1, 0.01)
ys <- dnorm(xs, 0, 0.4 / 2.57)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples2[['phi_sck[1]']], 50,
flim = c(0,1),
display_name=expression(phi[schock]))
util$plot_expectand_pushforward(samples2[['phi_sck[2]']], 50,
flim = c(0,1),
display_name=expression(phi[schock]),
add = T)
xs <- seq(0, 1, 0.01)
ys <- dbeta(xs, 2, 20)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples2[['omega_conc_sck']], 50,
flim = c(0.5,1),
display_name=expression(omega[shock]^concordant))
xs <- seq(0.5, 1, 0.01)
ys <- dbeta(xs, 230, 14)
lines(xs, ys, lwd=2.5, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
util$plot_expectand_pushforward(samples2[['omega_nonconc_sck']], 50,
flim = c(0,0.5),
display_name=expression(omega[shock]^nonconcordant))
xs <- seq(0, 1, 0.01)
ys <- dbeta(xs, 1, 20)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
So, what parameters in this fit (orange) are different than the previous fit with only one stand (green)?
\tau_\text{outer}^\text{shock} seems a bit smaller, which would suggest it should better capture intermediate shocks, no..?
The probability of non-concordant tree shock \omega_\text{shock}^\text{non-concordant} (i.e. given the stand is not shocking), is also quite different. To have a different perspective, we can look at the two product of probabilities that lead to a shock. (Note: stand 2 is only available in the second fit)
In the second fit, I kind of feel that the difference between \phi_\text{shock} * \omega_\text{shock}^\text{concordant} (Site 1, Q50 = 0.16) and (1-\phi_\text{shock}) * \omega_\text{shock}^\text{non-concordant} (Site 1, Q50 = 0.06) is not as high as I would expect. In other words, at the tree level, the proportion of non-concordant shocks relative to concordant shocks is higher than I would expect.
12 A step back to stay on track
I probably rushed things and skipped a step! I implemented the censoring only after the Cauchy implementation, when it should have been the very first thing to do. Boo me. Below, results with a normally distributed shock model WITH censoring.
data <- readRDS(file = file.path(wd, 'output', 'model', 'data_15nov2025_az592_az628.rds'))
data$rw_obs <- ifelse(data$log_rw_obs == log(1e-4), 0, exp(data$log_rw_obs))
if(FALSE){
fit <- stan(file=file.path(wd, 'model', 'stan', 'model9_marginalshocks_zerotauinner_censoring.stan'),
data=data, seed=5838293, chains = 4,
warmup=1000, iter=2024, refresh=10)
saveRDS(fit, file.path(wd, 'model/shocks/output/marginal', 'fit_normal_censoring.rds'))
}else{
fit <- readRDS(file.path(wd, 'model/shocks/output/marginal', 'fit_normal_censoring.rds'))
}Woa! A flood of tree depth warnings.
diagnostics <- util$extract_hmc_diagnostics(fit)
util$check_all_hmc_diagnostics(diagnostics) Chain 1: 1024 of 1024 transitions (100.00%)
saturated the maximum treedepth of 10.
Chain 2: 1024 of 1024 transitions (100.00%)
saturated the maximum treedepth of 10.
Chain 3: 1024 of 1024 transitions (100.00%)
saturated the maximum treedepth of 10.
Chain 4: 1024 of 1024 transitions (100.00%)
saturated the maximum treedepth of 10.
Numerical trajectories that saturate the maximum treedepth have
terminated prematurely. Increasing max_depth above 10 should result in
more expensive, but more efficient, Hamiltonian transitions.samples <- util$extract_expectand_vals(fit)
base_samples <- util$filter_expectands(samples,
c('alpha',
'beta_gdd', 'beta_sm', 'beta_vpd',
'rho_sh', 'gamma_sh',
'rho_sp', 'gamma_sp',
'sigma',
'phi_sck', 'omega_conc_sck'),
check_arrays = TRUE)
util$check_all_expectand_diagnostics(base_samples)rho_sp[1]:
Chain 1: hat{ESS} (71.699) is smaller than desired (100).
Chain 2: hat{ESS} (72.291) is smaller than desired (100).
Chain 3: hat{ESS} (93.643) is smaller than desired (100).
Chain 4: hat{ESS} (80.025) is smaller than desired (100).
Small empirical effective sample sizes result in imprecise Markov chain
Monte Carlo estimators.Let’s still do some retrodictive checks! On the first stand:
par(mfcol = c(4,2))
for(t in 1:4){
idxs_t <- data$tree_start_idxs[t]:data$tree_end_idxs[t]
names <- sapply(idxs_t,
function(x) paste0('log_rw_pred[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='log_rw_pred', main = paste0('Tree ', data$uniq_tree_ids[t]),
display_ylim=c(-8, 2), display_xlim = c(1980, 2010))
points(data$years[idxs_t], ifelse(data$log_rw_obs[idxs_t] == log(1e-4), NA, data$log_rw_obs[idxs_t]), pch=16, cex=1, col="white")
points(data$years[idxs_t], ifelse(data$log_rw_obs[idxs_t] == log(1e-4), NA, data$log_rw_obs[idxs_t]), pch=16, cex=0.6, col="black")
abline(v = c(2002,1996), lty = 3, col="darkgrey")
names <- sapply(idxs_t,
function(x) paste0('sck_state[', x, ']'))
ticks <- c(unlist(lapply(seq(1980, 2000,10), function(x) c(x, rep(NA, 9)))), '2010')
util$plot_disc_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab='Shock state',
display_ylim=c(-0.1, 1.1),
xticklabs = ticks)
abline(v = c(2002,1996), lty = 3, col="darkgrey")
names <- sapply(idxs_t,
function(x) paste0('delta_sck[', x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$years[idxs_t],
xlab="", ylab=expression(delta[shock]),
display_ylim=c(-8, 3), display_xlim = c(1980, 2010))
abline(v = c(2002,1996), lty = 3, col="darkgrey")
names <- sapply(1:data$N_all_years,
function(x) paste0('f_sh[', data$stand_idxs[t],',',x, ']'))
util$plot_conn_pushforward_quantiles(samples, names, data$all_years,
xlab="", ylab='Stand GP',
display_ylim=c(-8, 3), display_xlim = c(1980, 2010))
abline(v = c(2002,1996), lty = 3, col="darkgrey")
}
The model seems to over-regularize the large negative shocks. We obtain a \tau_\text{outer}^\text{shock} rather small compared to what our broad prior would suggest:
par(mfcol = c(1,2))
util$plot_expectand_pushforward(samples[['outer_tau_sck']], 50,
flim = c(0,5),
display_name=expression(tau[shock]^outer))
plot(1, type = "n", xlab = "", ylab = "", xlim = c(0, 10), ylim = c(0, 0.2),
xaxs = "i", yaxs = "i", bty = "n", yaxt = "n")
xs <- seq(0, 10, 0.01)
ys <- dnorm(xs, 0, 10 / 2.57)
lines(xs, ys, lwd=2.7, col='white')
lines(xs, ys, lwd=1.7, col=util$c_mid_teal)
This would lead me to consider a distribution a bit less rigid.